Micron Document




Algorithme de Chudnovski
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
top
L' algorithme de Chudnovsky est une méthode rapide pour calculer les chiffres de π, basée sur les formules de π de Ramanujan. Il a été publié par les frères Chudnovsky en 1988.

Il a été utilisé pour le calcul du record mondial de 2 700 milliards de chiffres de π en décembre 2009, 10 000 milliards de chiffres en octobre 2011, 22 400 milliards de chiffres en novembre 2016cite-ref-1[1], 31 400 milliards de chiffres en septembre 2018. – janvier 2019cite-ref-2[2], 50 000 milliards de chiffres le 29 janvier 2020cite-ref-3[3], 62 800 milliards de chiffres le 14 août 2021cite-ref-4[4], 100 000 milliards de chiffres le 21 mars 2022cite-ref-5[5], et 105 000 milliards de chiffres le 14 mars 2024cite-ref-6[6].

Contents


──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Algorithme

L'algorithme est basé sur l'opposé du nombre de Heegner d = − − 163 {\displaystyle d=-163} , la fonction j j ( 1 + i 163 2 ) = − − 640320 3 {\displaystyle j\left({\tfrac {1+i{\sqrt {163}}}{2}}\right)=-640320^{3}} , et sur la série hypergéométrique généralisée à convergence rapide ci-dessous: 1 π π = 12 ∑ ∑ k = 0 ∞ ∞ ( − − 1 ) k ( 6 k ) ! ( 545140134 k + 13591409 ) ( 3 k ) ! ( k ! ) 3 ( 640320 ) 3 k + 3 / 2 {\displaystyle {\frac {1}{\pi }}=12\sum _{k=0}^{\infty }{\frac {(-1)^{k}(6k)!(545140134k+13591409)}{(3k)!(k!)^{3}(640320)^{3k+3/2}}}} Une preuve détaillée de cette formule peut être trouvée icicite-ref-7[7] :

Cette identité est similaire à certaines formules de Ramanujan impliquant π, et est un exemple de série Ramanujan-Sato.

La complexité temporelle de l'algorithme est en O ( n ( log ⁡ ⁡ n ) 3 ) {\displaystyle O\left(n(\log n)^{3}\right)} cite-ref-8[8].

Optimisations

La technique d'optimisation utilisée pour les calculs du record du monde est appelée division binaire.

Fractionnement binaire

Un facteur de 1 / 640320 3 / 2 {\textstyle 1/{640320^{3/2}}} peut être retiré de la somme et simplifié en 1 π π = 1 426880 10005 ∑ ∑ k = 0 ∞ ∞ ( − − 1 ) k ( 6 k ) ! ( 545140134 k + 13591409 ) ( 3 k ) ! ( k ! ) 3 ( 640320 ) 3 k {\displaystyle {\frac {1}{\pi }}={\frac {1}{426880{\sqrt {10005}}}}\sum _{k=0}^{\infty }{\frac {(-1)^{k}(6k)!(545140134k+13591409)}{(3k)!(k!)^{3}(640320)^{3k}}}}

Soit f ( n ) = ( − − 1 ) n ( 6 n ) ! ( 3 n ) ! ( n ! ) 3 ( 640320 ) 3 n {\displaystyle f(n)={\frac {(-1)^{n}(6n)!}{(3n)!(n!)^{3}(640320)^{3n}}}} , et en substituant dans la somme. 1 π π = 1 426880 10005 ∑ ∑ k = 0 ∞ ∞ f ( k ) ⋅ ⋅ ( 545140134 k + 13591409 ) {\displaystyle {\frac {1}{\pi }}={\frac {1}{426880{\sqrt {10005}}}}\sum _{k=0}^{\infty }{f(k)\cdot (545140134k+13591409)}} f ( n ) f ( n − − 1 ) {\displaystyle {\frac {f(n)}{f(n-1)}}} peut être simplifié en − − ( 6 n − − 1 ) ( 2 n − − 1 ) ( 6 n − − 5 ) 10939058860032000 n 3 {\displaystyle {\frac {-(6n-1)(2n-1)(6n-5)}{10939058860032000n^{3}}}} , donc f ( n ) = f ( n − − 1 ) ⋅ ⋅ − − ( 6 n − − 1 ) ( 2 n − − 1 ) ( 6 n − − 5 ) 10939058860032000 n 3 {\displaystyle f(n)=f(n-1)\cdot {\frac {-(6n-1)(2n-1)(6n-5)}{10939058860032000n^{3}}}} f ( 0 ) = 1 {\displaystyle f(0)=1} par définition de f {\displaystyle f} , donc f ( n ) = ∏ ∏ j = 1 n − − ( 6 j − − 1 ) ( 2 j − − 1 ) ( 6 j − − 5 ) 10939058860032000 j 3 {\displaystyle f(n)=\prod _{j=1}^{n}{\frac {-(6j-1)(2j-1)(6j-5)}{10939058860032000j^{3}}}} Cette définition de f {\displaystyle f} n'est pas défini pour n = 0 {\displaystyle n=0} , on calcule donc le premier terme de la somme et l'ajoute à la nouvelle définition de f {\displaystyle f} en partant de n = 1 {\displaystyle n=1} 1 π π = 1 426880 10005 ( 13591409 + ∑ ∑ k = 1 ∞ ∞ ( ∏ ∏ j = 1 k − − ( 6 j − − 1 ) ( 2 j − − 1 ) ( 6 j − − 5 ) 10939058860032000 j 3 ) ⋅ ⋅ ( 545140134 k + 13591409 ) ) {\displaystyle {\frac {1}{\pi }}={\frac {1}{426880{\sqrt {10005}}}}{\Bigg (}13591409+\sum _{k=1}^{\infty }{{\Bigg (}\prod _{j=1}^{k}{\frac {-(6j-1)(2j-1)(6j-5)}{10939058860032000j^{3}}}{\Bigg )}\cdot (545140134k+13591409)}{\Bigg )}} Soit P ( a , b ) = ∏ ∏ j = a b − − 1 − − ( 6 j − − 1 ) ( 2 j − − 1 ) ( 6 j − − 5 ) {\displaystyle P(a,b)=\prod _{j=a}^{b-1}{-(6j-1)(2j-1)(6j-5)}} et Q ( a , b ) = ∏ ∏ j = a b − − 1 10939058860032000 j 3 {\displaystyle Q(a,b)=\prod _{j=a}^{b-1}{10939058860032000j^{3}}} , donc 1 π π = 1 426880 10005 ( 13591409 + ∑ ∑ k = 1 ∞ ∞ P ( 1 , k + 1 ) Q ( 1 , k + 1 ) ⋅ ⋅ ( 545140134 k + 13591409 ) ) {\displaystyle {\frac {1}{\pi }}={\frac {1}{426880{\sqrt {10005}}}}{\Bigg (}13591409+\sum _{k=1}^{\infty }{{\frac {P(1,k+1)}{Q(1,k+1)}}\cdot (545140134k+13591409)}{\Bigg )}} Soit S ( a , b ) = ∑ ∑ k = a b − − 1 P ( a , k + 1 ) Q ( a , k + 1 ) ⋅ ⋅ ( 545140134 k + 13591409 ) {\displaystyle S(a,b)=\sum _{k=a}^{b-1}{{\frac {P(a,k+1)}{Q(a,k+1)}}\cdot (545140134k+13591409)}} et R ( a , b ) = Q ( a , b ) ⋅ ⋅ S ( a , b ) {\displaystyle R(a,b)=Q(a,b)\cdot S(a,b)} π π = 426880 10005 13591409 + S ( 1 , ∞ ∞ ) {\displaystyle \pi ={\frac {426880{\sqrt {10005}}}{13591409+S(1,\infty )}}} S ( 1 , ∞ ∞ ) {\displaystyle S(1,\infty )} est une limite qui ne peut être qu'approché, on calcule à la place S ( 1 , n ) {\displaystyle S(1,n)} et lorsque n {\displaystyle n} se rapproches de ∞ ∞ {\displaystyle \infty } , l'approximation de π π {\displaystyle \pi } l'approximation s'améliore. π π ≈ ≈ 426880 10005 13591409 + S ( 1 , n ) {\displaystyle \pi \approx {\frac {426880{\sqrt {10005}}}{13591409+S(1,n)}}} Par définition originale de R {\displaystyle R} , S ( a , b ) = R ( a , b ) Q ( a , b ) {\displaystyle S(a,b)={\frac {R(a,b)}{Q(a,b)}}} π π ≈ ≈ 426880 10005 ⋅ ⋅ Q ( 1 , n ) 13591409 Q ( 1 , n ) + R ( 1 , n ) {\displaystyle \pi \approx {\frac {426880{\sqrt {10005}}\cdot Q(1,n)}{13591409Q(1,n)+R(1,n)}}}

Calcul récursif des fonctions

• P ( a , b ) = P ( a , m ) ⋅ ⋅ P ( m , b ) {\displaystyle P(a,b)=P(a,m)\cdot P(m,b)}

• Q ( a , b ) = Q ( a , m ) ⋅ ⋅ Q ( m , b ) {\displaystyle Q(a,b)=Q(a,m)\cdot Q(m,b)}
• S ( a , b ) = S ( a , m ) + P ( a , m ) Q ( a , m ) S ( m , b ) {\displaystyle S(a,b)=S(a,m)+{\frac {P(a,m)}{Q(a,m)}}S(m,b)}
• R ( a , b ) = Q ( m , b ) R ( a , m ) + P ( a , m ) R ( m , b ) {\displaystyle R(a,b)=Q(m,b)R(a,m)+P(a,m)R(m,b)}

Construction de la récursion

Si b = a + 1 {\displaystyle b=a+1}

• P ( a , a + 1 ) = − − ( 6 a − − 1 ) ( 2 a − − 1 ) ( 6 a − − 5 ) {\displaystyle P(a,a+1)=-(6a-1)(2a-1)(6a-5)}
• Q ( a , a + 1 ) = 10939058860032000 a 3 {\displaystyle Q(a,a+1)=10939058860032000a^{3}}
• S ( a , a + 1 ) = P ( a , a + 1 ) Q ( a , a + 1 ) ⋅ ⋅ ( 545140134 a + 13591409 ) {\displaystyle S(a,a+1)={\frac {P(a,a+1)}{Q(a,a+1)}}\cdot (545140134a+13591409)}
• R ( a , a + 1 ) = P ( a , a + 1 ) ⋅ ⋅ ( 545140134 a + 13591409 ) {\displaystyle R(a,a+1)=P(a,a+1)\cdot (545140134a+13591409)}

Code Python

import decimal
def binary_split(a, b):
if b == a + 1:
Pab = -(6*a - 5)*(2*a - 1)*(6*a - 1)
Qab = 10939058860032000 * a**3
Rab = Pab * (545140134*a + 13591409)
else:
m = (a + b) // 2
Pam, Qam, Ram = binary_split(a, m)
Pmb, Qmb, Rmb = binary_split(m, b)
Pab = Pam * Pmb
Qab = Qam * Qmb
Rab = Qmb * Ram + Pam * Rmb
return Pab, Qab, Rab
def chudnovsky(n):
"""Chudnovsky algorithm."""
P1n, Q1n, R1n = binary_split(1, n)
return (426880 * decimal.Decimal(10005).sqrt() * Q1n) / (13591409*Q1n + R1n)
print(chudnovsky(2)) # 3.141592653589793238462643384
decimal.getcontext().prec = 100
for n in range(2,10):
print(f"{n=} {chudnovsky(n)}") # 3.14159265358979323846264338...

Remarques

e π π 163 ≈ ≈ 640320 3 + 743.99999999999925 … … {\displaystyle e^{\pi {\sqrt {163}}}\approx 640320^{3}+743.99999999999925\dots }
640320 3 / 24 = 10939058860032000 {\displaystyle 640320^{3}/24=10939058860032000}
545140134 = 163 ⋅ ⋅ 127 ⋅ ⋅ 19 ⋅ ⋅ 11 ⋅ ⋅ 7 ⋅ ⋅ 3 2 ⋅ ⋅ 2 {\displaystyle 545140134=163\cdot 127\cdot 19\cdot 11\cdot 7\cdot 3^{2}\cdot 2}
13591409 = 13 ⋅ ⋅ 1045493 {\displaystyle 13591409=13\cdot 1045493}

Voir aussi

• Série Ramanujan-Sato
• L'algorithme de Borwein

Références

• (en) Cet article est partiellement ou en totalité issu de l’article de Wikipédia en anglais intitulé « Chudnovsky algorithm » (voir la liste des auteurs).

cite-note-11. « 22.4 Trillion Digits of Pi », www.numberworld.org
cite-note-22. « Google Cloud Topples the Pi Record », www.numberworld.org/
cite-note-33. « The Pi Record Returns to the Personal Computer », www.numberworld.org/
cite-note-44. « Pi-Challenge - Weltrekordversuch der FH Graubünden - FH Graubünden », www.fhgr.ch (consulté le 17 août 2021)
cite-note-55. « Calculating 100 trillion digits of pi on Google Cloud », cloud.google.com (consulté le 10 juin 2022)
cite-note-66. yee2024Yee, « Limping to a new Pi Record of 105 Trillion Digits », NumberWorld.org, 14 mars 2024 (consulté le 16 mars 2024)
cite-note-77. milla2018lorenz-milla2018Lorenz Milla, « A detailed proof of the Chudnovsky formula with means of basic complex analysis », 2018 (arXiv 1809.00533)
cite-note-88. « y-cruncher - Formulas », www.numberworld.org (consulté le 25 février 2018)

• Portail de l'analyse